Introducción

Si es su primera vez trabajando con los datos, no está muy familiarizado con el paquete o simplemente quiere conocer más sobre el proyecto y el objetivo de estos ejemplos y el paquete verdata, consulte: https://github.com/HRDAG/CO-examples/blob/main/Introducción/output/Introducción.html antes de continuar.

En este ejemplo se ilustrará el proceso de estimación del total de víctimas por año del hecho (1990–2017).

Importando resultados estratificación

En este ejemplo se evidenció el proceso para la estratificación en la que agrupamos víctimas con características similares. Teniendo este insumo procedimos a guardar los resultados en nuestra maquina local o computador. En este ejemplo retomaremos dicho input para continuar con el siguiente paso que es -precisamente- la estimación:

results <- readRDS(here::here("Resultados-CEV/output-estratificacion/reclutamiento-yy_hecho-edad_c.rds"))

Estimación víctimas por año del hecho

Ahora, definidos los estratos, se calculan las estimaciones para esta violación en particular; para este paso se usa la función mse del paquete verdata. Esta función permite preparar los datos, revisar si hay estimaciones precalculadas para ese estrato y estimar, en caso de que no las haya. Esta función tomará como insumo la información de las fuentes, es decir, aquellas columnas que comienzan por in_. Para que un estrato sea estimable, se requiere un mínimo de 3 fuentes válidas. Si un estrato no es estimable la función arrojará NA.

Como se mencionó, considerando que el proceso de estimación toma tiempo y recursos computacionales, esta función le permite usar estimaciones ya calculadas en el proyecto. Para esto, usted debe descargar las estimaciones publicadas acá.

mse(
  stratum_data,
  stratum_name,
  estimates_dir = NULL,
  min_n = 1,
  K = NULL,
  buffer_size = 10000,
  sampler_thinning = 1000,
  seed = 19481210,
  burnin = 10000,
  n_samples = 10000,
  posterior_thinning = 500
)

Algunos de los argumentos de esta función se explican de la siguiente forma¹

• stratum_data: Data frame que incluye la información del estrato de interés (data frames que creamos antes).

• stratum_name: Es el nombre del estrato.

• estimates_dir: Es la ruta (opcional) o el path de la carpeta o archivo que contiene las estimaciones precalculadas. Esto le permite a la función buscar entre las estimacines si el estrato que usted quiere analizar ya fue estimado.

estimaciones_dir <- here::here("estimaciones")

Al final el resultado será un data frame con cinco columnas: una columna que indica si el estrato es válido o no (TRUE o FALSE); el número de muestras N de la distribución posterior (NA si el estrato no es válido); las fuentes válidas que se usaron en la estimación valid_sources; el número de observaciones de las listas que son válidas del estrato (n_obs) y el nombre del estrato stratum_name. Si un estrato es estimable, este data frame va a contener 1000 muestras para cada estrato.

Antes de esto es importante aclarar que este proceso se hará a través de una iteración (repetir el proceso), ya que este proceso se aplicará a cada una de nuestras réplicas. Es decir:

En primera instancia guardaremos nuestros resultados por réplica en un objeto llamado estimates. Posteriormente procederemos a usar nuestra iteración en cada réplica, siendo réplicas el objeto definido como replicas <- paste0("R", 1:10). Ahora, map2_dfr nos será bastante útil ya que permite aplicar nuestra función de mse a cada combinación de nuestras estratificaciones. Es decir, esta función permite calcular la función de mse a cada una de nuestras combinaciones de datos: strata_data y las variables de estratificación stratification_vars. Por último agregaremos la columna de replica para ubicar de qué versión o réplica de los datos viene nuestra información y con do.call pegaremos nuestras bases, es decir, las réplicas. Esto lo veremos a continuación:

estimates <- list() 

replicas <- paste0("R", 1:10)

for (replica in replicas) {
    
    strata_data <- results[[replica]]$strata_data
    stratification_vars <- results[[replica]]$stratification_vars
    
    estimate_mse <- purrr::map2_dfr(.x = strata_data,
                                    .y = stratification_vars,
                                    .f = ~ {
                                      mse <- verdata::mse(stratum_data = .x,
                                                          stratum_name = .y,
                                                          estimates_dir = estimaciones_dir)
                               mse$replica <- replica
                               mse
                           })
    
    estimates[[replica]] <- estimate_mse
}

union_mse <- do.call(rbind, estimates)

union_mse <- union_mse %>%
    mutate(stratum_name = paste(pull(stratum_name, 1),
                                pull(stratum_name, 2),
                                sep = "-"))

paged_table(union_mse, options = list(rows.print = 10, cols.print = 8))

Estos son los resultados de mse. Vemos las cinco (5) primeras columnas mencionadas en el que cabe destacar que todos los estratos son válidos, es decir, no vemos ningún NA en nuestra tabla; lo que indica que en cada estrato hay al menos 3 listas con al menos 1 víctima. Además, vemos nuestro N el cual evidencia 1000 muestras aleatorias de la distribución posterior de la cantidad de víctimas probables para cada estrato y réplica.

Combinación de las estimaciones

Ahora, ya que tenemos esta información, nuestro último paso en el trabajo de estimación para víctimas -desagregadas por año del hecho- es combinar nuestras estimaciones. Esto lo haremos a través de nuestra función combine_estimates el cual permite obtener los intervalos creibles (no de confianza). Pero antes, debemos realizar una pequeña transformación a nuestra anterior tabla union_mse, esto con el fin de obtener nuestros resultados desagregados esta característica (y proceder luego a combinar). Esto lo presentaremos a continuación:

tabla_sampler <- union_mse %>%
    separate(stratum_name, into = c("yy_hecho", "edad_c"), 
             sep = "-", remove = FALSE, extra = "merge") %>% 
    rename(replicate = replica) %>% 
    group_by(stratum_name, replicate) %>% 
    mutate(sample_number = glue("sample_{row_number()}")) %>% 
    pivot_wider(id_cols = c("replicate", "stratum_name", "yy_hecho", "n_obs"),
                values_from = N,
                names_from = sample_number) 

paged_table(tabla_sampler, options = list(rows.print = 10, cols.print = 5))

En primera instancia vemos que la función separate del paquete tidyr nos separa la columna de stratum_name en las variables originales de perpetrador y año del hecho, luego agrupamos por nuestra columna stratum_name y cada una de nuestras 10 réplicas.

Adicionalmente generamos una columna denominada sample_number en el que {row_number()} corresponde a cada fila de N. Es decir, articulando esto con nuestra siguiente línea de código, vemos que pivot_wider convierte nuestro dataframe de un formato largo a ancho, o, en otras palabras, los valores que están en N se “trasladan” a las columnas que comienzan por sample_. Como habíamos explicado anteriormente, cada estrato para la réplica (R1) contiene 1000 muestras aleatorias de la distribución posterior, por lo que cada valor de este N se expanderá a cada columna de sample_ y por tanto, esas serán 1000 columnas, desde sample_1 hasta sample_1000.

Teniendo esta transformación, procederemos a agrupar por nuestra variable de interés (año del hecho) y réplica, seguido de la suma de estas muestras de la distribución posterior de cada estrato que forma parte de este agrupamiento.

tabla_agrupacion <- tabla_sampler %>% 
    group_by(yy_hecho, replicate) %>% 
    summarize(across(c(starts_with("sample_"), "n_obs"), sum)) 

paged_table(tabla_agrupacion, options = list(rows.print = 10, cols.print = 5))

En otras palabras, primero agrupamos por departamento del hecho y réplica, es decir -por ejemplo- agrupar por el año de 1990 y la réplica 1, 1990 y réplica 10, y así sucesivamente. Cuando tengamos este tipo de agrupación (como lo vemos en la tabla_agrupacion) tomaremos -por ejemplo- la réplica 1 (“R1”) y la categoría o año de 1990 y sumaremos los valores para el sample_ (es como tomar: ejemplo <- tabla_sampler %>% filter(replicate == "R1", yy_hecho == "1990") y sumar: sum(ejemplo$sample_1), dándonos lo que nos está arrojando la primera fila de la columna sample_1, es decir, 196).

Sumado a esto, volveremos a transformar nuestros datos de ancho a largo, es decir, volveremos a una estructura muy parecida a tabla_sampler, pero esta vez nuestro N muestra la suma de las muestras de la distribución posterior por año del hecho y réplica que hicimos en el paso anterior:

tiempo_agrupacion <- tabla_agrupacion %>%
    group_by(yy_hecho) %>%
    pivot_longer(starts_with("sample_"), 
                 names_to = "replicate_num",
                 values_to = "N") %>%
    ungroup() %>%
    select(-replicate_num) %>%
    rename(stratum_name = yy_hecho) 

paged_table(tiempo_agrupacion, options = list(rows.print = 10, cols.print = 5))

Para finalizar, vemos que esta información está desagregada por año del hecho y réplica, pero, como sabemos, deseamos ver los resultados únicamente por año del hecho. Por tal razón agruparemos esta información por dicha variable (que ahora se denomina stratum_name):

final_agrupacion <- tiempo_agrupacion %>%
    group_by(stratum_name)

estimates_tabla <- final_agrupacion %>%
    group_split() %>%
    map_dfr(.f = combine_estimates) %>%
    bind_cols(group_keys(final_agrupacion)) %>% 
    select(stratum_name, N_025, N_mean, N_975) %>% 
    rename(yy_hecho = stratum_name)

paged_table(estimates_tabla, options = list(rows.print = 10, cols.print = 5))

Vemos pues que, después de agrupar, tenemos los resultados derivados de estimates_tabla cuyo código permite separar nuestras categorías (o o años), luego utilizamos la función de map_dfr en el que podemos aplicar combine_estimates a cada una de nuestras categorías. Por último -con bind_cols- unimos estos resultados con las variables de agrupación, es decir, con yy_hecho que la denominamos stratum_name.

Seguido de esto importaremos nuestros resultados a nivel documentado e imputado:

tabla_doc_imp <- arrow::read_parquet(here::here("Resultados-CEV/output-imputados/reclutamiento-anio-imputados.parquet"))

paged_table(tabla_doc_imp, options = list(rows.print = 10, cols.print = 5))

Por último uniremos esta base con nuestros resultados combinados y graficaremos nuestros resultados:

estimates_final <- estimates_tabla %>% 
    mutate(yy_hecho = as.character(yy_hecho))

tabla_final <- dplyr::left_join(estimates_tabla, tabla_doc_imp)

tabla_final <- tabla_final %>% 
    arrange(desc(N_mean))
    
paged_table(tabla_final, options = list(rows.print = 10, cols.print = 5))

n_warn <- 2000

combine_estimates_year <- tabla_final %>%
    mutate(max_var = case_when(
      (N_975 > n_warn) ~ n_warn,
      TRUE ~ NA_real_)) %>%
    mutate(N_975 = ifelse(!is.na(max_var), max_var, N_975))

combine_estimates_year <- combine_estimates_year %>% 
    arrange(desc(N_mean)) %>% 
    mutate(yy_hecho = as.numeric(yy_hecho))

mr_observed_ttl <- glue("Observado")
mr_replicates_ttl <- glue("Imputado")
mr_universo_ttl <- glue("Estimado")

g <- combine_estimates_year %>%
  ggplot(aes(x = yy_hecho)) +
  geom_line(aes(y = observed,
                fill = mr_observed_ttl), color = "black", size = 1) +
  geom_line(aes(y = imp_mean,  fill = mr_replicates_ttl),color = "#1F74B1", 
            show.legend = FALSE, size = 1) +
  geom_ribbon(aes(ymin = N_025, ymax = N_975, fill = mr_universo_ttl),
              alpha = 0.5) +
  geom_point(aes(y = max_var), pch = 21, color = 'darkgreen', fill = "green",
             size = 1, stroke = 2) +
  theme_minimal() +
  xlab("") +
  ylab("Número de víctimas") +
  theme(axis.text.x = element_text(size = 11, angle = 90),
        axis.title.y = element_text(size = 11),
        axis.ticks.x = element_line(size = 0.1)) +
  scale_x_continuous(breaks = seq(1985, 2020, 5)) +
  theme(legend.position = "bottom") +
  scale_fill_manual(values = c("darkgreen", "#1F74B1", "black" ), name = "")

print(g)

El proyecto JEP-CEV-HRDAG estimó que para el año 2 013 el número de víctimas estuvo entre 1 669 y 4 781 con una probabilidad del 95% (intervalo de credibilidad). En otras palabras, hay una alta probabilidad (95% de credibilidad) de que el número de de victimas en este año se encuentre dentro de este rango, mostrando a su vez el valor más probable de 2 561 victimas. Adicionalmente se puede evidenciar que la varianza es tan alta que no se presenta en la gráfica.